%     FASTENER DOCTOR - программа для лечения и конвертации крепежа
%     Copyright (C) 2018  Михаил Козлов
%     Сайт проекта: <http://www.fastenerdoctor.ru/>
%     Адрес электронной почты для связи с разработчиком:
%     <doctor.fastener@gmail.com>
% 
%     Эта программа является свободным программным обеспечением: Вы можете
%     распространять ее и/или изменять, соблюдая условия Универсальной общественной
%     лицензии GNU, опубликованной Фондом свободного программного обеспечения,
%     либо редакции 3, либо (на Ваше усмотрение) любой редакции, выпущенной позже.
% 
%     Эта программа распространяется в надежде на то, что она окажется полезной,
%     но БЕЗ КАКИХ-ЛИБО ГАРАНТИЙ, включая даже подразумеваемую гарантию ТОВАРНОЙ ПРИГОДНОСТИ,
%     либо ПРИГОДНОСТИ ДЛЯ ОПРЕДЕЛЕННЫХ ЦЕЛЕЙ. Ознакомьтесь с Универсальной общественной
%     лицензией GNU для получения более подробной информации.
% 
%     Вы должны были получить копию Универсальной общественной лицензии GNU
%     вместе с этой программой. Если Вы ее не получили, то перейдите по адресу:
%     <http://www.gnu.org/licenses/>.
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.

function [A1,B1,D1, A2,B2,D2] = Laminate_Stiffness_Matrices(N_Plies, Mat_Prop_Table, Thickness_Table, Stacking_Angle_Table)

% Формировать ли матрицы жесткости ламината относительно нейтральной базовой плоскости

Account_For_Neutral_Plane = 'TRUE'; % TRUE - в качестве базовой плоскости по каждому направлению будет использована нейтральная базовая плоскость вдоль этого направления
                                    % FALSE - в качестве базовой плоскости по каждому направлению будет использована срединная плоскость ламината

% Массивы синусов и косинусов

s = sind(Stacking_Angle_Table);
c = cosd(Stacking_Angle_Table);

% Поворот жесткостных свойств слоев

Q = Mat_Prop_Table;
Q_lam = zeros(N_Plies,6);

for i = 1:N_Plies
    
    Q_lam(i,1) = Q(i,1)*c(i)^4 + Q(i,4)*s(i)^4 + 2*(Q(i,2) + 2*Q(i,6))*s(i)^2*c(i)^2 - 4*s(i)*c(i)*(c(i)^2*Q(i,3) + s(i)^2*Q(i,5));
    Q_lam(i,2) = (Q(i,1) + Q(i,4) - 4*Q(i,6))*s(i)^2*c(i)^2 + Q(i,2)*(c(i)^4 + s(i)^4) + 2*s(i)*c(i)*(Q(i,3) - Q(i,5))*(c(i)^2 - s(i)^2);
    Q_lam(i,3) = (Q(i,1) - Q(i,2) - 2*Q(i,6))*s(i)*c(i)^3 + (Q(i,2) - Q(i,4) + 2*Q(i,6))*s(i)^3*c(i) + (c(i)^2 - s(i)^2)*(c(i)^2*Q(i,3) + s(i)^2*Q(i,5)) + 2*s(i)^2*c(i)^2*(Q(i,5) - Q(i,3));
    Q_lam(i,4) = Q(i,1)*s(i)^4 + Q(i,4)*c(i)^4 + 2*(Q(i,2) + 2*Q(i,6))*s(i)^2*c(i)^2 + 4*s(i)*c(i)*(s(i)^2*Q(i,3) + c(i)^2*Q(i,5));
    Q_lam(i,5) = (Q(i,1) - Q(i,2) - 2*Q(i,6))*s(i)^3*c(i) + (Q(i,2) - Q(i,4) + 2*Q(i,6))*s(i)*c(i)^3 + (c(i)^2 - s(i)^2)*(s(i)^2*Q(i,3) + c(i)^2*Q(i,5)) + 2*s(i)^2*c(i)^2*(Q(i,3) - Q(i,5));
    Q_lam(i,6) = (Q(i,1) + Q(i,4) - 2*Q(i,2) - 2*Q(i,6))*s(i)^2*c(i)^2 + Q(i,6)*(s(i)^4 + c(i)^4) + 2*s(i)*c(i)*(c(i)^2 - s(i)^2)*(Q(i,3) - Q(i,5));
    
end

% Массив t_k

t_k = Thickness_Table;

% Расчет координат центров слоев

ksi = zeros(N_Plies,1);

ksi(1) = t_k(1)/2;

for i=2:N_Plies
    ksi(i) = ksi(i-1) + (t_k(i) + t_k(i-1))/2;
end

% Расчет положения нейтральной плоскости при смятии вдоль оси 1

temp_1 = 0;
temp_2 = 0;

for i=1:N_Plies
    temp_1 = temp_1 + t_k(i)*Q_lam(i,1)*ksi(i);
    temp_2 = temp_2 + t_k(i)*Q_lam(i,1);
end

c1 = temp_1/temp_2;

% Заполнение массивов p_k_1, q_k_1

z_k_1 = zeros(N_Plies+1,1);

if (strcmp(Account_For_Neutral_Plane, 'TRUE'))
    z_k_1(1) = -c1;
else
    z_k_1(1) = -sum(t_k)/2;
end

for i = 2:N_Plies+1
    z_k_1(i) = z_k_1(i-1) + t_k(i-1);
end

p_k_1 = zeros(N_Plies,1);
q_k_1 = zeros(N_Plies,1);

for i = 1:N_Plies
    
    p_k_1(i) = z_k_1(i+1)^2 - z_k_1(i)^2;
    q_k_1(i) = z_k_1(i+1)^3 - z_k_1(i)^3;
    
end

% Сборка матрицы жесткости ламината при смятии вдоль оси 1

A1 = zeros(3,3);
B1 = zeros(3,3);
D1 = zeros(3,3);

for i = 1:N_Plies
    
    A1(1,1) = A1(1,1) + Q_lam(i,1) * t_k(i);
    A1(1,2) = A1(1,2) + Q_lam(i,2) * t_k(i);
    A1(1,3) = A1(1,3) + Q_lam(i,3) * t_k(i);
    A1(2,2) = A1(2,2) + Q_lam(i,4) * t_k(i);
    A1(2,3) = A1(2,3) + Q_lam(i,5) * t_k(i);
    A1(3,3) = A1(3,3) + Q_lam(i,6) * t_k(i);
    
    B1(1,1) = B1(1,1) + Q_lam(i,1) * p_k_1(i);
    B1(1,2) = B1(1,2) + Q_lam(i,2) * p_k_1(i);
    B1(1,3) = B1(1,3) + Q_lam(i,3) * p_k_1(i);
    B1(2,2) = B1(2,2) + Q_lam(i,4) * p_k_1(i);
    B1(2,3) = B1(2,3) + Q_lam(i,5) * p_k_1(i);
    B1(3,3) = B1(3,3) + Q_lam(i,6) * p_k_1(i);
    
    D1(1,1) = D1(1,1) + Q_lam(i,1) * q_k_1(i);
    D1(1,2) = D1(1,2) + Q_lam(i,2) * q_k_1(i);
    D1(1,3) = D1(1,3) + Q_lam(i,3) * q_k_1(i);
    D1(2,2) = D1(2,2) + Q_lam(i,4) * q_k_1(i);
    D1(2,3) = D1(2,3) + Q_lam(i,5) * q_k_1(i);
    D1(3,3) = D1(3,3) + Q_lam(i,6) * q_k_1(i);
    
end

B1 = B1/2;
D1 = D1/3;

A1(2,1) = A1(1,2);
A1(3,1) = A1(1,3);
A1(3,2) = A1(2,3);

B1(2,1) = B1(1,2);
B1(3,1) = B1(1,3);
B1(3,2) = B1(2,3);

D1(2,1) = D1(1,2);
D1(3,1) = D1(1,3);
D1(3,2) = D1(2,3);

% Расчет положения нейтральной плоскости при смятии вдоль оси 2

temp_1 = 0;
temp_2 = 0;

for i=1:N_Plies
    temp_1 = temp_1 + t_k(i)*Q_lam(i,4)*ksi(i);
    temp_2 = temp_2 + t_k(i)*Q_lam(i,4);
end

c2 = temp_1/temp_2;

% Заполнение массивов p_k_2, q_k_2

z_k_2 = zeros(N_Plies+1,1);

if (strcmp(Account_For_Neutral_Plane, 'TRUE'))
    z_k_2(1) = -c2;
else
    z_k_2(1) = -sum(t_k)/2;
end

for i = 2:N_Plies+1
    z_k_2(i) = z_k_2(i-1) + t_k(i-1);
end

p_k_2 = zeros(N_Plies,1);
q_k_2 = zeros(N_Plies,1);

for i = 1:N_Plies
    
    p_k_2(i) = z_k_2(i+1)^2 - z_k_2(i)^2;
    q_k_2(i) = z_k_2(i+1)^3 - z_k_2(i)^3;
    
end

% Сборка матрицы жесткости ламината при смятии вдоль оси 2

A2 = zeros(3,3);
B2 = zeros(3,3);
D2 = zeros(3,3);

for i = 1:N_Plies
    
    A2(1,1) = A2(1,1) + Q_lam(i,1) * t_k(i);
    A2(1,2) = A2(1,2) + Q_lam(i,2) * t_k(i);
    A2(1,3) = A2(1,3) + Q_lam(i,3) * t_k(i);
    A2(2,2) = A2(2,2) + Q_lam(i,4) * t_k(i);
    A2(2,3) = A2(2,3) + Q_lam(i,5) * t_k(i);
    A2(3,3) = A2(3,3) + Q_lam(i,6) * t_k(i);
    
    B2(1,1) = B2(1,1) + Q_lam(i,1) * p_k_2(i);
    B2(1,2) = B2(1,2) + Q_lam(i,2) * p_k_2(i);
    B2(1,3) = B2(1,3) + Q_lam(i,3) * p_k_2(i);
    B2(2,2) = B2(2,2) + Q_lam(i,4) * p_k_2(i);
    B2(2,3) = B2(2,3) + Q_lam(i,5) * p_k_2(i);
    B2(3,3) = B2(3,3) + Q_lam(i,6) * p_k_2(i);
    
    D2(1,1) = D2(1,1) + Q_lam(i,1) * q_k_2(i);
    D2(1,2) = D2(1,2) + Q_lam(i,2) * q_k_2(i);
    D2(1,3) = D2(1,3) + Q_lam(i,3) * q_k_2(i);
    D2(2,2) = D2(2,2) + Q_lam(i,4) * q_k_2(i);
    D2(2,3) = D2(2,3) + Q_lam(i,5) * q_k_2(i);
    D2(3,3) = D2(3,3) + Q_lam(i,6) * q_k_2(i);
    
end

B2 = B2/2;
D2 = D2/3;

A2(2,1) = A2(1,2);
A2(3,1) = A2(1,3);
A2(3,2) = A2(2,3);

B2(2,1) = B2(1,2);
B2(3,1) = B2(1,3);
B2(3,2) = B2(2,3);

D2(2,1) = D2(1,2);
D2(3,1) = D2(1,3);
D2(3,2) = D2(2,3);
